home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume2 / aplictns / matlab / doc.2 < prev    next >
Internet Message Format  |  1988-11-03  |  56KB
  1. Path: xanth!nic.MR.NET!tank!uwvax!rutgers!mailrus!ulowell!page
  2. From: page@swan.ulowell.edu (Bob Page)
  3. Newsgroups: comp.sources.amiga
  4. Subject: v02i051:  matlab - matrix laboratory, Part11/11 (doc02/02)
  5. Message-ID: <10026@swan.ulowell.edu>
  6. Date: 2 Nov 88 21:55:51 GMT
  7. Organization: University of Lowell, Computer Science Dept.
  8. Lines: 1746
  9. Approved: page@swan.ulowell.edu
  10.  
  11. Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
  12. Posting-number: Volume 2, Issue 51
  13. Archive-name: applications/matlab/doc.2
  14.  
  15. #    This is a shell archive.
  16. #    Remove everything above and including the cut line.
  17. #    Then run the rest of the file through sh.
  18. #----cut here-----cut here-----cut here-----cut here----#
  19. #!/bin/sh
  20. # shar:    Shell Archiver
  21. #    Run the following text with /bin/sh to create:
  22. #    doc.2
  23. # This archive created: Wed Nov  2 16:14:44 1988
  24. cat << \SHAR_EOF > doc.2
  25.    HESS(A)         - same as EIG.
  26.  
  27.  
  28.  
  29.  
  30.  
  31.  
  32.  
  33.  
  34.  
  35.  
  36.  
  37. MATLAB, page 35
  38.  
  39.  
  40.  
  41.      Minor modifications were made to all these subroutines.  The
  42. LINPACK  routines  were  changed  to  replace the Fortran complex
  43. arithmetic with explicit references to real and imaginary  parts.
  44. Since  most of the floating point arithmetic is concentrated in a
  45. few low-level subroutines which perform  vector  operations  (the
  46. Basic  Linear  Algebra  Subprograms),  this  was not an extensive
  47. change.  It also facilitated implementation of the FLOP and  CHOP
  48. features  which count and optionally truncate each floating point
  49. operation.
  50.  
  51.      The EISPACK subroutine COMQR2 was modified to  allow  access
  52. to  the  Schur  triangular  form, ordinarily just an intermediate
  53. result.   IMTQL2  was  modified  to  make  computation   of   the
  54. eigenvectors   optional.    Both  subroutines  were  modified  to
  55. eliminate the machine-dependent accuracy parameter  and  all  the
  56. EISPACK subroutines were changed to include FLOP and CHOP.
  57.  
  58.      The algorithms employed for the  POLY  and  ROOTS  functions
  59. illustrate  an  interesting  aspect  of  the  modern  approach to
  60. eigenvalue computation.   POLY(A)  generates  the  characteristic
  61. polynomial  of  A  and  ROOTS(POLY(A))  finds  the  roots of that
  62. polynomial, which are, of course, the eigenvalues of A . But both
  63. POLY  and  ROOTS  use  EISPACK eigenvalues subroutines, which are
  64. based on similarity transformations.  So the  classical  approach
  65. which  characterizes  eigenvalues  as roots of the characteristic
  66. polynomial is actually reversed.
  67.  
  68.      If A is an n by n matrix, POLY(A) produces the  coefficients
  69. C(1) through C(n+1), with C(1) = 1, in
  70.  
  71.       DET(z*EYE-A) = C(1)*z**n + ... + C(n)*z + C(n+1) .
  72.  
  73. The algorithm can be expressed compactly using MATLAB:
  74.  
  75.       Z = EIG(A);
  76.       C = 0*ONES(n+1,1);  C(1) = 1;
  77.       for j = 1:n, C(2:j+1) = C(2:j+1) - Z(j)*C(1:j);
  78.       C
  79.  
  80. This recursion is easily derived by expanding the product
  81.  
  82.       (z - z(1))*(z - z(2))* ... * (z-z(n)) .
  83.  
  84. It is possible to prove that POLY(A) produces the coefficients in
  85. the  characteristic  polynomial of a matrix within roundoff error
  86. of  A .  This is true even if the  eigenvalues  of  A  are  badly
  87. conditioned.    The  traditional  algorithms  for  obtaining  the
  88. characteristic polynomial which do not use the eigenvalues do not
  89. have such satisfactory numerical properties.
  90.  
  91.      If C is a vector with n+1  components,  ROOTS(C)  finds  the
  92. roots of the polynomial of degree n ,
  93.  
  94.  
  95.  
  96.  
  97.  
  98.  
  99.  
  100.  
  101.  
  102.  
  103. MATLAB, page 36
  104.  
  105.  
  106.  
  107.        p(z) = C(1)*z**n + ... + C(n)*z + C(n+1) .
  108.  
  109. The algorithm simply involves computing the  eigenvalues  of  the
  110. companion matrix:
  111.  
  112.       A = 0*ONES(n,n)
  113.       for j = 1:n, A(1,j) = -C(j+1)/C(1);
  114.       for i = 2:n, A(i,i-1) = 1;
  115.       EIG(A)
  116.  
  117. It is possible to prove that the results produced are  the  exact
  118. eigenvalues  of  a  matrix within roundoff error of the companion
  119. matrix A, but this does not mean that they are the exact roots of
  120. a  polynomial with coefficients within roundoff error of those in
  121. C .  There are more accurate, more efficient methods for  finding
  122. polynomial  roots,  but  this  approach has the crucial advantage
  123. that it does not require very much additional code.
  124.  
  125.      The elementary functions EXP, LOG, SQRT, SIN, COS  and  ATAN
  126. are  applied  to  square  matrices  by  diagonalizing the matrix,
  127. applying the functions to the  individual  eigenvalues  and  then
  128. transforming back.  For example, EXP(A) is computed by
  129.  
  130.       <X,D> = EIG(A);
  131.       for j = 1:n, D(j,j) = EXP(D(j,j));
  132.       X*D/X
  133.  
  134. This is essentially method number 14  out  of  the  19  'dubious'
  135. possibilities described in [8].  It is dubious because it doesn't
  136. always work.  The matrix of eigenvectors  X  can  be  arbitrarily
  137. badly  conditioned  and  all  accuracy lost in the computation of
  138. X*D/X.  A warning message is printed if RCOND(X) is  very  small,
  139. but  this  only  catches the extreme cases.  An example of a case
  140. not detected is when A has a double eigenvalue, but theoretically
  141. only  one  linearly  independent  eigenvector associated with it.
  142. The computed eigenvalues will be separated by  something  on  the
  143. order  of the square root of the roundoff level.  This separation
  144. will be reflected in RCOND(X) which will probably  not  be  small
  145. enough to trigger the error message.  The computed EXP(A) will be
  146. accurate to only half precision.  Better methods  are  known  for
  147. computing EXP(A), but they do not easily extend to the other five
  148. functions and would require a considerable amount  of  additional
  149. code.
  150.  
  151.      The expression A**p is evaluated by repeated  multiplication
  152. if p is an integer greater than 1.  Otherwise it is evaluated by
  153.  
  154.       <X,D> = EIG(A);
  155.       for j = 1:n, D(j,j) = EXP(p*LOG(D(j,j)))
  156.       X*D/X
  157.  
  158. This suffers from the same potential loss of  accuracy  if  X  is
  159. badly conditioned.  It was partly for this reason that the case p
  160.  
  161.  
  162.  
  163.  
  164.  
  165.  
  166.  
  167.  
  168.  
  169. MATLAB, page 37
  170.  
  171.  
  172.  
  173. = 1 is included in the general case.  Comparison of A**1  with  A
  174. gives some idea of the loss of accuracy for other values of p and
  175. for the elementary functions.
  176.  
  177.      RREF, the reduced row echelon form, is of some  interest  in
  178. theoretical  linear algebra, although it has little computational
  179. value.  It is included in MATLAB for  pedagogical  reasons.   The
  180. algorithm  is essentially Gauss-Jordan elimination with detection
  181. of negligible columns applied to rectangular matrices.
  182.  
  183.      There are three separate places in MATLAB where the rank  of
  184. a  matrix  is  implicitly  computed:  in RREF(A), in A\B for non-
  185. square A, and in  the  pseudoinverse  PINV(A).   Three  different
  186. algorithms  with  three  different criteria for negligibility are
  187. used and so it is possible that three different values  could  be
  188. produced for the same matrix.  With RREF(A), the rank of A is the
  189. number of nonzero rows.  The elimination algorithm used for  RREF
  190. is  the  fastest of the three rank-determining algorithms, but it
  191. is the least sophisticated numerically and  the  least  reliable.
  192. With  A\B,  the  algorithm  is  essentially  that used by example
  193. subroutine SQRST  in  chapter  9  of  the  LINPACK  guide.   With
  194. PINV(A),   the   algorithm   is   based  on  the  singular  value
  195. decomposition and is described  in  chapter  11  of  the  LINPACK
  196. guide.   The  SVD  algorithm  is the most time-consuming, but the
  197. most reliable and is therefore also used for RANK(A).
  198.  
  199.      The  uniformly  distributed  random  numbers  in  RAND   are
  200. obtained  from  the  machine-independent  random number generator
  201. URAND described in [9].  It is possible  to  switch  to  normally
  202. distributed   random   numbers,   which   are  obtained  using  a
  203. transformation also described in [9].
  204.  
  205.      The computation of
  206.  
  207.                 2    2
  208.           sqrt(a  + b )
  209.  
  210. is  required  in  many  matrix  algorithms,  particularly   those
  211. involving  complex  arithmetic.   A  new approach to carrying out
  212. this operation is described by Moler and Morrison [10].  It is  a
  213. cubically  convergent  algorithm  which  starts with  a  and  b ,
  214. rather than with their squares, and  thereby  avoids  destructive
  215. arithmetic underflows and overflows.  In MATLAB, the algorithm is
  216. used for complex modulus, Euclidean vector norm, plane rotations,
  217. and  the  shift  calculation in the eigenvalue and singular value
  218. iterations.
  219.  
  220.  
  221. 12.  FLOP and CHOP
  222.  
  223.      Detailed information about the amount of  work  involved  in
  224. matrix  calculations  and  the  resulting accuracy is provided by
  225. FLOP and CHOP.  The basic unit of work is the "flop", or floating
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235. MATLAB, page 38
  236.  
  237.  
  238.  
  239. point operation.  Roughly, one flop is one execution of a Fortran
  240. statement like
  241.  
  242.       S = S + X(I)*Y(I)
  243.  
  244. or
  245.  
  246.       Y(I) = Y(I) + T*X(I)
  247.  
  248. In other words, it consists of one floating point multiplication,
  249. together  with  one  floating  point  addition and the associated
  250. indexing and storage reference operations.
  251.  
  252.      MATLAB will  print  the  number  of  flops  required  for  a
  253. particular statement when the statement is terminated by an extra
  254. comma.  For example, the line
  255.  
  256.       n = 20;  RAND(n)*RAND(n);,
  257.  
  258. ends with an extra comma.  Two  20  by  20  random  matrices  are
  259. generated  and  multiplied  together.   The result is assigned to
  260. ANS, but the semicolon suppresses its printing.  The only  output
  261. is
  262.  
  263.         8800 flops
  264.  
  265. This is  n**3 + 2*n**2  flops,  n**2  for each random matrix  and
  266. n**3 for the product.
  267.  
  268.      FLOP is a predefined vector with two components.  FLOP(1) is
  269. the number of flops used by the most recently executed statement,
  270. except that statements with zero flops are ignored.  For example,
  271. after executing the previous statement,
  272.  
  273.       flop(1)/n**3
  274.  
  275. results in
  276.  
  277.       ANS   =
  278.  
  279.           1.1000
  280.  
  281.  
  282.      FLOP(2) is the cumulative total of all the flops used  since
  283. the beginning of the MATLAB session.  The statement
  284.  
  285.       FLOP = <0 0>
  286.  
  287. resets the total.
  288.  
  289.      There are several difficulties  associated  with  keeping  a
  290. precise  count  of  floating  point  operations.  An  addition or
  291. subtraction that is not paired with a multiplication  is  usually
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301. MATLAB, page 39
  302.  
  303.  
  304.  
  305. counted as a flop. The same is true of an isolated multiplication
  306. that is  not  paired  with  an  addition.   Each  floating  point
  307. division counts as a flop.  But the number of operations required
  308. by system dependent library functions such as square root  cannot
  309. be  counted, so most elementary functions are arbitrarily counted
  310. as using only one flop.
  311.  
  312.      The  biggest  difficulty  occurs  with  complex  arithmetic.
  313. Almost  all operations on the real parts of matrices are counted.
  314. However, the operations on the  complex  parts  of  matrices  are
  315. counted only when they involve nonzero elements.  This means that
  316. simple operations on nonreal matrices require only about twice as
  317. many  flops as the same operations on real matrices.  This factor
  318. of two is not necessarily an accurate  measure  of  the  relative
  319. costs of real and complex arithmetic.
  320.  
  321.      The result of each floating  point  operation  may  also  be
  322. "chopped" to simulate a computer with a shorter word length.  The
  323. details of this chopping operation depend upon the format of  the
  324. floating point word.  Usually, the fraction in the floating point
  325. word  can  be  regarded  as  consisting  of  several   octal   or
  326. hexadecimal digits.  The least significant of these digits can be
  327. set to zero by a logical masking operation.  Thus the statement
  328.  
  329.       CHOP(p)
  330.  
  331. causes the  p  least significant octal or hexadecimal  digits  in
  332. the  result  of  each floating point operation to be set to zero.
  333. For example, if the computer being  used  has  an  IBM  360  long
  334. floating  point  word with 14 hexadecimal digits in the fraction,
  335. then CHOP(8) results in simulation of  a  computer  with  only  6
  336. hexadecimal  digits  in the fraction, i.e. a short floating point
  337. word. On a computer such as the CDC 6600 with  16  octal  digits,
  338. CHOP(8)  results in about the same accuracy because the remaining
  339. 8 octal digits represent the same number of bits as 6 hexadecimal
  340. digits.
  341.  
  342.      Some idea of the effect of CHOP on any particular system can
  343. be obtained by executing the following statements.
  344.  
  345.       long,   t = 1/10
  346.       long z, t = 1/10
  347.       chop(8)
  348.       long,   t = 1/10
  349.       long z, t = 1/10
  350.  
  351.  
  352.      The following Fortran subprograms illustrate more details of
  353. FLOP  and CHOP. The first subprogram is a simplified example of a
  354. system-dependent function used within MATLAB itself.  The  common
  355. variable  FLP  is essentially the first component of the variable
  356. FLOP.  The common variable CHP is initially zero, but it  is  set
  357. to  p  by the statement  CHOP(p).  To shorten the DATA statement,
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367. MATLAB, page 40
  368.  
  369.  
  370.  
  371. we assume there are only 6 hexadecimal digits.  We also assume an
  372. extension  of  Fortran  that  allows .AND. to be used as a binary
  373. operation between two real variables.
  374.  
  375.       REAL FUNCTION FLOP(X)
  376.       REAL X
  377.       INTEGER FLP,CHP
  378.       COMMON FLP,CHP
  379.       REAL MASK(5)
  380.       DATA MASK/ZFFFFFFF0,ZFFFFFF00,ZFFFFF000,ZFFFF0000,ZFFF00000/
  381.       FLP = FLP + 1
  382.       IF (CHP .EQ. 0) FLOP = X
  383.       IF (CHP .GE. 1 .AND. CHP .LE. 5) FLOP = X .AND. MASK(CHP)
  384.       IF (CHP .GE. 6) FLOP = 0.0
  385.       RETURN
  386.       END
  387.  
  388.  
  389.      The following subroutine illustrates a typical  use  of  the
  390. previous  function  within MATLAB.  It is a simplified version of
  391. the Basic Linear Algebra Subprogram that adds a  scalar  multiple
  392. of  one  vector  to another.  We assume here that the vectors are
  393. stored with a memory increment of one.
  394.  
  395.       SUBROUTINE SAXPY(N,TR,TI,XR,XI,YR,YI)
  396.       REAL TR,TI,XR(N),XI(N),YR(N),YI(N),FLOP
  397.       IF (N .LE. 0) RETURN
  398.       IF (TR .EQ. 0.0 .AND. TI .EQ. 0.0) RETURN
  399.       DO 10 I = 1, N
  400.          YR(I) = FLOP(YR(I) + TR*XR(I) - TI*XI(I))
  401.          YI(I) = YI(I) + TR*XI(I) + TI*XR(I)
  402.          IF (YI(I) .NE. 0.0D0) YI(I) = FLOP(YI(I))
  403.    10 CONTINUE
  404.       RETURN
  405.       END
  406.  
  407.  
  408.      The  saxpy  operation  is  perhaps  the   most   fundamental
  409. operation  within  LINPACK.  It is used in the computation of the
  410. LU, the QR and the  SVD  factorizations,  and  in  several  other
  411. places.   We  see  that  adding  a  multiple of one vector with n
  412. components to another uses n flops if the vectors  are  real  and
  413. between  n  and  2*n  flops if the vectors have nonzero imaginary
  414. components.
  415.  
  416.      The permanent MATLAB variable EPS is reset by the  statement
  417. CHOP(p).   Its new value is usually the smallest inverse power of
  418. two that satisfies the Fortran logical test
  419.  
  420.             FLOP(1.0+EPS) .GT. 1.0
  421.  
  422. However, if EPS had been directly reset to a  larger  value,  the
  423. old value is retained.
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433. MATLAB, page 41
  434.  
  435.  
  436.  
  437.  
  438.  
  439. 13.  Communicating with other programs
  440.  
  441.      There  are  four  different  ways  MATLAB  can  be  used  in
  442. conjunction with other programs:
  443.       -- USER,
  444.       -- EXEC,
  445.       -- SAVE and LOAD,
  446.       -- MATZ, CALL and RETURN .
  447.  
  448.      Let us illustrate each of  these  by  the  following  simple
  449. example.
  450.  
  451.       n = 6
  452.       for i = 1:n, for j = 1:n, a(i,j) = abs(i-j);
  453.       A
  454.       X = inv(A)
  455.  
  456.  
  457.      The example  A  could be introduced into MATLAB  by  writing
  458. the following Fortran subroutine.
  459.  
  460.          SUBROUTINE USER(A,M,N,S,T)
  461.          DOUBLE PRECISION A(1),S,T
  462.          N = IDINT(A(1))
  463.          M = N
  464.          DO 10 J = 1, N
  465.          DO 10 I = 1, N
  466.             K = I + (J-1)*M
  467.             A(K) = IABS(I-J)
  468.       10 CONTINUE
  469.          RETURN
  470.          END
  471.  
  472. This subroutine should be compiled  and  linked  into  MATLAB  in
  473. place   of  the  original  version  of  USER.   Then  the  MATLAB
  474. statements
  475.  
  476.       n = 6
  477.       A = user(n)
  478.       X = inv(A)
  479.  
  480. do the job.
  481.  
  482.      The example A could be generated by  storing  the  following
  483. text in a file named, say, EXAMPLE .
  484.  
  485.       for i = 1:n, for j = 1:n, a(i,j) = abs(i-j);
  486.  
  487. Then the MATLAB statements
  488.  
  489.       n = 6
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  
  499. MATLAB, page 42
  500.  
  501.  
  502.  
  503.       exec('EXAMPLE',0)
  504.       X = inv(A)
  505.  
  506. have the desired effect.  The 0 as the optional second  parameter
  507. of exec indicates that the text in the file should not be printed
  508. on the terminal.
  509.  
  510.      The matrices A and X could also be  stored  in  files.   Two
  511. separate main programs would be involved.  The first is:
  512.  
  513.          PROGRAM MAINA
  514.          DOUBLE PRECISION A(10,10)
  515.          N = 6
  516.          DO 10 J = 1, N
  517.          DO 10 I = 1, N
  518.             A(I,J) = IABS(I-J)
  519.       10 CONTINUE
  520.          OPEN(UNIT=1,FILE='A')
  521.          WRITE(1,101) N,N
  522.      101 FORMAT('A   ',2I4)
  523.          DO 20 J = 1, N
  524.             WRITE(1,102) (A(I,J),I=1,N)
  525.       20 CONTINUE
  526.      102 FORMAT(4Z18)
  527.          END
  528.  
  529. The OPEN statement may take different forms on different systems.
  530. It  attaches  Fortran  logical unit number 1 to the file named A.
  531. The FORMAT  number  102  may  also  be  system  dependent.   This
  532. particular one is appropriate for hexadecimal computers with an 8
  533. byte double precision floating point  word.   Check,  or  modify,
  534. MATLAB subroutine SAVLOD.
  535.  
  536.      After this program is executed, enter MATLAB  and  give  the
  537. following statements:
  538.  
  539.       load('A')
  540.       X = inv(A)
  541.       save('X',X)
  542.  
  543. If all goes according to plan, this will read the matrix  A  from
  544. the  file A, invert it, store the inverse in X and then write the
  545. matrix X on the file X .  The following program can then access X
  546. .
  547.  
  548.          PROGRAM MAINX
  549.          DOUBLE PRECISION X(10,10)
  550.          OPEN(UNIT=1,FILE='X')
  551.          REWIND 1
  552.          READ (1,101) ID,M,N
  553.      101 FORMAT(A4,2I4)
  554.          DO 10 J = 1, N
  555.             READ(1,102) (X(I,J),I=1,M)
  556.  
  557.  
  558.  
  559.  
  560.  
  561.  
  562.  
  563.  
  564.  
  565. MATLAB, page 43
  566.  
  567.  
  568.  
  569.       10 CONTINUE
  570.      102 FORMAT(4Z18)
  571.          ...
  572.          ...
  573.  
  574.  
  575.      The most elaborate mechanism  involves  using  MATLAB  as  a
  576. subroutine within another program.  Communication with the MATLAB
  577. stack is accomplished using subroutine MATZ which is  distributed
  578. with  MATLAB,  but  which  is  not  used  by  MATLAB itself.  The
  579. preample of MATZ is:
  580.  
  581.       SUBROUTINE MATZ(A,LDA,M,N,IDA,JOB,IERR)
  582.       INTEGER LDA,M,N,IDA(1),JOB,IERR
  583.       DOUBLE PRECISION A(LDA,N)
  584. C
  585. C     ACCESS MATLAB VARIABLE STACK
  586. C     A IS AN M BY N MATRIX, STORED IN AN ARRAY WITH
  587. C         LEADING DIMENSION LDA.
  588. C     IDA IS THE NAME OF A.
  589. C         IF IDA IS AN INTEGER K LESS THAN 10, THEN THE NAME IS 'A'K
  590. C         OTHERWISE, IDA(1:4) IS FOUR CHARACTERS, FORMAT 4A1.
  591. C     JOB =  0  GET REAL A FROM MATLAB,
  592. C         =  1  PUT REAL A INTO MATLAB,
  593. C         = 10  GET IMAG PART OF A FROM MATLAB,
  594. C         = 11  PUT IMAG PART OF A INTO MATLAB.
  595. C     RETURN WITH NONZERO IERR AFTER MATLAB ERROR MESSAGE.
  596. C
  597. C     USES MATLAB ROUTINES STACKG, STACKP AND ERROR
  598.  
  599.  
  600.      The preample of subroutine MATLAB is:
  601.  
  602.  
  603.       SUBROUTINE MATLAB(INIT)
  604. C     INIT = 0 FOR FIRST ENTRY, NONZERO FOR SUBSEQUENT ENTRIES
  605.  
  606.  
  607.      To do our example, write the following program:
  608.  
  609.          DOUBLE PRECISION A(10,10),X(10,10)
  610.          INTEGER IDA(4),IDX(4)
  611.          DATA LDA/10/
  612.          DATA IDA/'A',' ',' ',' '/, IDX/'X',' ',' ',' '/
  613.          CALL MATLAB(0)
  614.          N = 6
  615.          DO 10 J = 1, N
  616.          DO 10 I = 1, N
  617.             A(I,J) = IABS(I-J)
  618.       10 CONTINUE
  619.          CALL MATZ(A,LDA,N,N,IDA,1,IERR)
  620.          IF (IERR .NE. 0) GO TO ...
  621.          CALL MATLAB(1)
  622.  
  623.  
  624.  
  625.  
  626.  
  627.  
  628.  
  629.  
  630.  
  631. MATLAB, page 44
  632.  
  633.  
  634.  
  635.          CALL MATZ(X,LDA,N,N,IDX,0,IERR)
  636.          IF (IERR .NE. 0) GO TO ...
  637.          ...
  638.          ...
  639.  
  640. When this program is executed, the call to MATLAB(0) produces the
  641. MATLAB greeting, then waits for input.  The command
  642.  
  643.          return
  644.  
  645. sends control back to our  example  program.   The  matrix  A  is
  646. generated  by the program and sent to the stack by the first call
  647. to MATZ.  The call to MATLAB(1) produces the MATLAB prompt.  Then
  648. the statements
  649.  
  650.          X = inv(A)
  651.          return
  652.  
  653. will invert our matrix, put the result on the stack and  go  back
  654. to our program.  The second call to MATZ will retrieve X .
  655.  
  656.      By the way, this matrix  X  is interesting. Take a  look  at
  657. round(2*(n-1)*X).
  658.  
  659.  
  660.  
  661.  
  662. Acknowledgement.
  663.  
  664.  
  665.      Most of the work on MATLAB  has  been  carried  out  at  the
  666. University  of  New  Mexico,  where  it is being supported by the
  667. National Science Foundation. Additional work has been done during
  668. visits  to  Stanford  Linear Accelerator Center, Argonne National
  669. Laboratory and Los Alamos Scientific  Laboratory,  where  support
  670. has been provided by NSF and the Department of Energy.
  671.  
  672.  
  673. References
  674.  
  675. [1]  J. J. Dongarra, J. R. Bunch, C. B. Moler and G. W.  Stewart,
  676.      LINPACK  Users'  Guide,  Society  for Industrial and Applied
  677.      Mathematics, Philadelphia, 1979.
  678.  
  679. [2]  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S.  Garbow,  Y.
  680.      Ikebe, V. C. Klema, C. B. Moler, Matrix Eigensystem Routines
  681.      -- EISPACK Guide, Lecture Notes in Computer Science,  volume
  682.      6, second edition, Springer-Verlag, 1976.
  683.  
  684. [3]  B. S. Garbow, J. M. Boyle, J.  J.  Dongarra,  C.  B.  Moler,
  685.      Matrix  Eigensystem  Routines  --  EISPACK  Guide Extension,
  686.      Lecture Notes in  Computer  Science,  volume  51,  Springer-
  687.      Verlag, 1977.
  688.  
  689.  
  690.  
  691.  
  692.  
  693.  
  694.  
  695.  
  696.  
  697. MATLAB, page 45
  698.  
  699.  
  700.  
  701. [4]  S. Cohen and  S.  Piper,  SPEAKEASY  III  Reference  Manual,
  702.      Speakeasy Computing Corp., Chicago, Ill., 1979.
  703.  
  704. [5]  J. H. Wilkinson  and  C.  Reinsch,  Handbook  for  Automatic
  705.      Computation,  volume  II,  Linear  Algebra, Springer-Verlag,
  706.      1971.
  707.  
  708. [6]  Niklaus Wirth, Algorithms  +  Data  Structures  =  Programs,
  709.      Prentice-Hall, 1976.
  710.  
  711. [7]  H. B. Keller and D. Sachs, "Calculations of the Conductivity
  712.      of  a  Medium Containing Cylindrical Inclusions", J. Applied
  713.      Physics 35, 537-538, 1964.
  714.  
  715. [8]  C. B. Moler and C. F. Van Loan,  Nineteen  Dubious  Ways  to
  716.      Compute  the  Exponential  of a Matrix, SIAM Review 20, 801-
  717.      836, 1979.
  718.  
  719. [9]  G. E. Forsythe, M. A. Malcolm  and  C.  B.  Moler,  Computer
  720.      Methods for Mathematical Computations, Prentice-Hall, 1977.
  721.  
  722. [10] C. B. Moler and D. R. Morrison, "Replacing square  roots  by
  723.      Pythagorean   sums",  University  of  New  Mexico,  Computer
  724.      Science  Department,   technical   report,   submitted   for
  725.      publication, 1980.
  726.  
  727.  
  728.  
  729.  
  730.  
  731.  
  732.  
  733.  
  734.  
  735.  
  736.  
  737.  
  738.  
  739.  
  740.  
  741.  
  742.  
  743.  
  744.  
  745.  
  746.  
  747.  
  748.  
  749.  
  750.  
  751.  
  752.  
  753.  
  754.  
  755.  
  756.  
  757.  
  758.  
  759.  
  760.  
  761.  
  762.  
  763. MATLAB, page 46
  764.  
  765.  
  766.  
  767. Appendix.  The HELP document
  768.  
  769. NEWS  MATLAB NEWS dated May, 1981.
  770.       This describes recent or local changes.
  771.       The new features added since the November,  1980,  printing
  772.       of the Users' Guide include DIARY, EDIT, KRON, MACRO, PLOT,
  773.       RAT, TRIL, TRIU and six element-by-element operations:
  774.             .*   ./   .\   .*.   ./.   .\.
  775.       Some additional  capabilities  have  been  added  to  EXIT,
  776.       RANDOM, RCOND, SIZE and SVD.
  777.  
  778. INTRO Welcome to MATLAB.
  779.  
  780.       Here are a few sample statements:
  781.  
  782.       A = <1 2; 3 4>
  783.       b = <5 6>'
  784.       x = A\b
  785.       <V,D> = eig(A),  norm(A-V*D/V)
  786.       help \ , help eig
  787.       exec('demo',7)
  788.  
  789.       For more information, see the MATLAB Users' Guide which  is
  790.       contained in file ...  or may be obtained from ... .
  791.  
  792. HELP  HELP gives assistance.
  793.       HELP HELP obviously prints this message.
  794.       To see all the HELP messages, list the file ... .
  795.  
  796. <     < > Brackets used in forming vectors and matrices.
  797.       <6.9  9.64  SQRT(-1)>  is  a  vector  with  three  elements
  798.       separated  by  blanks.   <6.9,  9.64, sqrt(-1)> is the same
  799.       thing.  <1+I 2-I 3>  and  <1 +I 2 -I 3>  are not the  same.
  800.       The first has three elements, the second has five.
  801.       <11 12 13; 21 22 23>  is a 2 by 3 matrix .   The  semicolon
  802.       ends the first row.
  803.  
  804.       Vectors and matrices can be used inside < > brackets.
  805.       <A B; C>  is allowed if the number of rows  of   A   equals
  806.       the  number  of rows of  B  and the number of columns of  A
  807.       plus the number of columns of   B   equals  the  number  of
  808.       columns  of   C  .   This  rule  generalizes in a hopefully
  809.       obvious way to allow fairly complicated constructions.
  810.  
  811.       A = < >  stores an empty matrix in  A , thereby removing it
  812.       from the list of current variables.
  813.  
  814.       For the use of < and > on the left of  the  =  in  multiple
  815.       assignment statements, see LU, EIG, SVD and so on.
  816.  
  817.       In WHILE and IF clauses, <>  means  less  than  or  greater
  818.       than,  i.e.  not  equal, < means less than, > means greater
  819.       than, <= means less than or equal, >= means greater than or
  820.  
  821.  
  822.  
  823.  
  824.  
  825.  
  826.  
  827.  
  828.  
  829. MATLAB, page 47
  830.  
  831.  
  832.  
  833.       equal.
  834.  
  835.       For the use of > and < to delineate macros, see MACRO.
  836.  
  837. >     See < .  Also see MACRO.
  838.  
  839. (     ( ) Used to indicate precedence in  arithmetic  expressions
  840.       in  the  usual way.  Used to enclose arguments of functions
  841.       in the usual way.  Used to enclose  subscripts  of  vectors
  842.       and  matrices  in  a  manner somewhat more general than the
  843.       usual way.  If  X   and   V  are  vectors,  then   X(V)  is
  844.       <X(V(1)),  X(V(2)),  ...,  X(V(N))> .  The components of  V
  845.       are rounded to nearest integers and used as subscripts.  An
  846.       error  occurs  if  any  such  subscript  is  less than 1 or
  847.       greater than the dimension of  X .  Some examples:
  848.       X(3)  is the third element of  X .
  849.       X(<1 2 3>)  is the first three elements of  X .  So is
  850.       X(<SQRT(2), SQRT(3), 4*ATAN(1)>)  .
  851.       If  X  has  N  components,  X(N:-1:1) reverses them.
  852.       The same indirect subscripting is used in matrices.  If   V
  853.       has   M  components and  W  has  N  components, then A(V,W)
  854.       is the  M by N  matrix formed from the elements of A  whose
  855.       subscripts are the elements of  V  and  W .  For example...
  856.       A(<1,5>,:) = A(<5,1>,:)  interchanges rows 1 and 5 of  A .
  857.  
  858. )     See  ( .
  859.  
  860. =     Used in assignment statements and to mean equality in WHILE
  861.       and IF clauses.
  862.  
  863. .     Decimal point.  314/100, 3.14  and   .314E1   are  all  the
  864.       same.
  865.  
  866.       Element-by-element multiplicative operations  are  obtained
  867.       using  .*  ,  ./  , or .\ .  For example, C = A ./ B is the
  868.       matrix with elements  c(i,j) = a(i,j)/b(i,j) .
  869.  
  870.       Kronecker tensor products and quotients are  obtained  with
  871.       .*. , ./.  and .\. .  See KRON.
  872.  
  873.       Two or  more  points  at  the  end  of  the  line  indicate
  874.       continuation.    The   total  line  length  limit  is  1024
  875.       characters.
  876.  
  877. ,     Used to separate matrix subscripts and function  arguments.
  878.       Used  at  the  end  of  FOR, WHILE and IF clauses.  Used to
  879.       separate statements  in  multi-statement  lines.   In  this
  880.       situation,  it  may  be  replaced  by semicolon to suppress
  881.       printing.
  882.  
  883. ;     Used inside brackets to end rows.
  884.       Used after an expression or statement to suppress printing.
  885.       See SEMI.
  886.  
  887.  
  888.  
  889.  
  890.  
  891.  
  892.  
  893.  
  894.  
  895. MATLAB, page 48
  896.  
  897.  
  898.  
  899. \     Backslash or matrix left division.   A\B   is  roughly  the
  900.       same  as   INV(A)*B  , except it is computed in a different
  901.       way.  If  A  is an N by N matrix and  B  is a column vector
  902.       with  N  components, or a matrix with several such columns,
  903.       then X = A\B  is the solution to  the  equation   A*X  =  B
  904.       computed  by  Gaussian  elimination.   A warning message is
  905.       printed if  A is badly scaled or nearly singular.
  906.       A\EYE produces the inverse of  A .
  907.  
  908.       If  A  is an  M by N  matrix with  M < or > N  and  B  is a
  909.       column vector with  M  components, or a matrix with several
  910.       such columns, then  X = A\B  is the solution in  the  least
  911.       squares  sense  to  the under- or overdetermined system  of
  912.       equations A*X = B .  The  effective  rank,  K,  of   A   is
  913.       determined  from  the  QR  decomposition  with pivoting.  A
  914.       solution  X  is  computed  which  has  at  most  K  nonzero
  915.       components  per column.  If  K < N this will usually not be
  916.       the same solution as PINV(A)*B .
  917.       A\EYE produces a generalized inverse of  A .
  918.  
  919.       If A and B have the  same  dimensions,  then  A  .\  B  has
  920.       elements a(i,j)\b(i,j) .
  921.  
  922.       Also, see EDIT.
  923.  
  924. /     Slash or matrix right division.  B/A  is roughly  the  same
  925.       as  B*INV(A) .  More precisely,  B/A = (A'\B')' .  See \ .
  926.  
  927.       IF A and B have the  same  dimensions,  then  A  ./  B  has
  928.       elements a(i,j)/b(i,j) .
  929.  
  930.       Two or more slashes together on a line indicate  a  logical
  931.       end of line.  Any following text is ignored.
  932.  
  933. '     Transpose.  X'  is the complex conjugate transpose of  X  .
  934.       Quote.   'ANY  TEXT'   is a vector whose components are the
  935.       MATLAB internal codes for the characters.  A  quote  within
  936.       the text is indicated by two quotes.  See DISP and FILE .
  937.  
  938. +     Addition.  X + Y .  X and Y must have the same dimensions.
  939.  
  940. -     Subtraction.  X  -  Y  .   X  and  Y  must  have  the  same
  941.       dimensions.
  942.  
  943. *     Matrix multiplication, X*Y .  Any scalar (1  by  1  matrix)
  944.       may multiply anything.  Otherwise, the number of columns of
  945.       X must equal the number of rows of Y .
  946.  
  947.       Element-by-element multiplication is obtained with X .* Y .
  948.  
  949.       The Kronecker tensor product is denoted by X .*. Y .
  950.  
  951.       Powers.  X**p  is  X  to the   p   power.   p   must  be  a
  952.  
  953.  
  954.  
  955.  
  956.  
  957.  
  958.  
  959.  
  960.  
  961. MATLAB, page 49
  962.  
  963.  
  964.  
  965.       scalar.  If  X  is a matrix, see  FUN .
  966.  
  967. :     Colon.  Used in subscripts,  FOR  iterations  and  possibly
  968.       elsewhere.
  969.       J:K  is the same as  <J, J+1, ..., K>
  970.       J:K  is empty if  J > K .
  971.       J:I:K  is the same as  <J, J+I, J+2I, ..., K>
  972.       J:I:K  is empty if  I > 0 and J > K or if I < 0 and J < K .
  973.       The colon notation can be used to pick out  selected  rows,
  974.       columns and elements of vectors and matrices.
  975.       A(:)  is all the  elements  of  A,  regarded  as  a  single
  976.       column.
  977.       A(:,J)  is the  J-th  column of A
  978.       A(J:K)  is  A(J),A(J+1),...,A(K)
  979.       A(:,J:K)  is  A(:,J),A(:,J+1),...,A(:,K) and so on.
  980.       For the use of the colon in the FOR statement, See FOR .
  981.  
  982. ABS   ABS(X)  is the absolute value, or complex modulus,  of  the
  983.       elements of X .
  984.  
  985. ANS   Variable created automatically  when  expressions  are  not
  986.       assigned to anything else.
  987.  
  988. ATAN  ATAN(X)  is the arctangent of  X .  See FUN .
  989.  
  990. BASE  BASE(X,B) is a vector containing the base B  representation
  991.       of   X  .   This is often used in conjunction with DISPLAY.
  992.       DISPLAY(X,B)  is  the  same  as  DISPLAY(BASE(X,B)).    For
  993.       example,    DISP(4*ATAN(1),16)   prints   the   hexadecimal
  994.       representation of pi.
  995.  
  996. CHAR  CHAR(K)  requests  an  input  line  containing   a   single
  997.       character  to  replace  MATLAB  character  number  K in the
  998.       following table.  For example, CHAR(45) replaces backslash.
  999.       CHAR(-K) replaces the alternate character number K.
  1000.  
  1001.                 K  character alternate name
  1002.               0 - 9   0 - 9    0 - 9   digits
  1003.              10 - 35  A - Z    a - z   letters
  1004.                36                      blank
  1005.                37       (        (     lparen
  1006.                38       )        )     rparen
  1007.                39       ;        ;     semi
  1008.                40       :        |     colon
  1009.                41       +        +     plus
  1010.                42       -        -     minus
  1011.                43       *        *     star
  1012.                44       /        /     slash
  1013.                45       \        $     backslash
  1014.                46       =        =     equal
  1015.                47       .        .     dot
  1016.                48       ,        ,     comma
  1017.                49       '        "     quote
  1018.  
  1019.  
  1020.  
  1021.  
  1022.  
  1023.  
  1024.  
  1025.  
  1026.  
  1027. MATLAB, page 50
  1028.  
  1029.  
  1030.  
  1031.                50       <        [     less
  1032.                51       >        ]     great
  1033.  
  1034. CHOL  Cholesky factorization.  CHOL(X)  uses  only  the  diagonal
  1035.       and upper triangle of  X .  The lower triangular is assumed
  1036.       to be the (complex conjugate) transpose of the  upper.   If
  1037.       X   is  positive  definite,  then  R = CHOL(X)  produces an
  1038.       upper triangular  R  so that  R'*R = X .   If   X   is  not
  1039.       positive definite, an error message is printed.
  1040.  
  1041. CHOP  Truncate arithmetic.  CHOP(P) causes P places to be chopped
  1042.       off   after   each   arithmetic   operation  in  subsequent
  1043.       computations.  This means  P  hexadecimal  digits  on  some
  1044.       computers  and  P octal digits on others.  CHOP(0) restores
  1045.       full precision.
  1046.  
  1047. CLEAR Erases all variables, except EPS, FLOP, EYE and RAND.
  1048.       X = <>  erases only variable  X .  So does CLEAR X .
  1049.  
  1050. COND  Condition number in 2-norm.  COND(X) is the  ratio  of  the
  1051.       largest singular value of  X  to the smallest.
  1052.  
  1053. CONJG CONJG(X)  is the complex conjugate of  X .
  1054.  
  1055. COS   COS(X)  is the cosine of  X .  See FUN .
  1056.  
  1057. DET   DET(X)  is the determinant of the square matrix  X .
  1058.  
  1059. DIAG  If  V  is  a  row  or  column  vector  with  N  components,
  1060.       DIAG(V,K)   is a square matrix of order  N+ABS(K)  with the
  1061.       elements of  V  on the K-th diagonal.  K = 0  is  the  main
  1062.       diagonal,  K  >  0  is above the main diagonal and K < 0 is
  1063.       below the main diagonal.  DIAG(V)  simply puts  V   on  the
  1064.       main diagonal.
  1065.       eg. DIAG(-M:M) + DIAG(ONES(2*M,1),1) + DIAG(ONES(2*M,1),-1)
  1066.       produces a tridiagonal matrix of order 2*M+1 .
  1067.       IF  X  is a matrix,  DIAG(X,K)  is a column  vector  formed
  1068.       from the elements of the K-th diagonal of  X .
  1069.       DIAG(X)  is the main diagonal of  X .
  1070.       DIAG(DIAG(X))  is a diagonal matrix .
  1071.  
  1072. DIARY DIARY('file') causes a  copy  of  all  subsequent  terminal
  1073.       input and most of the resulting output to be written on the
  1074.       file. DIARY(0) turns it off.  See FILE.
  1075.  
  1076. DISP  DISPLAY(X) prints X  in  a  compact  format.   If  all  the
  1077.       elements  of  X  are  integers  between 0 and 51, then X is
  1078.       interpreted  as  MATLAB  text  and   printed   accordingly.
  1079.       Otherwise,  +  ,  -   and  blank  are printed for positive,
  1080.       negative and zero elements.  Imaginary parts are ignored.
  1081.       DISP(X,B) is the same as DISP(BASE(X,B)).
  1082.  
  1083. EDIT  There  are  no   editing   features   available   on   most
  1084.  
  1085.  
  1086.  
  1087.  
  1088.  
  1089.  
  1090.  
  1091.  
  1092.  
  1093. MATLAB, page 51
  1094.  
  1095.  
  1096.  
  1097.       installations and EDIT is not a command.  However, on a few
  1098.       systems a command line consisting of a single  backslash  \
  1099.       will  cause  the local file editor to be called with a copy
  1100.       of the  previous  input  line.   When  the  editor  returns
  1101.       control to MATLAB, it will execute the line again.
  1102.  
  1103. EIG   Eigenvalues and eigenvectors.
  1104.       EIG(X) is a vector containing the eigenvalues of  a  square
  1105.       matrix  X .
  1106.       <V,D>  =  EIG(X)   produces  a  diagonal  matrix    D    of
  1107.       eigenvalues  and  a  full  matrix  V  whose columns are the
  1108.       corresponding eigenvectors so that  X*V = V*D .
  1109.  
  1110. ELSE  Used with IF .
  1111.  
  1112. END   Terminates the scope  of  FOR,  WHILE  and  IF  statements.
  1113.       Without  END's,  FOR  and WHILE repeat all statements up to
  1114.       the end of the line.  Each END is paired with  the  closest
  1115.       previous  unpaired FOR or WHILE and serves to terminate its
  1116.       scope.  The line
  1117.       FOR I=1:N, FOR J=1:N, A(I,J)=1/(I+J-1); A
  1118.       would cause A to be printed  N**2  times, once for each new
  1119.       element.  On the other hand, the line
  1120.       FOR I=1:N, FOR J=1:N, A(I,J)=1/(I+J-1); END, END, A
  1121.       will lead to only the final printing of  A .
  1122.       Similar considerations apply to WHILE.
  1123.       EXIT terminates execution of loops or of MATLAB itself.
  1124.  
  1125. EPS   Floating point relative  accuracy.   A  permanent  variable
  1126.       whose  value is initially the distance from 1.0 to the next
  1127.       largest floating point number.  The  value  is  changed  by
  1128.       CHOP,  and  other values may be assigned.  EPS is used as a
  1129.       default tolerance by PINV and RANK.
  1130.  
  1131. EXEC  EXEC('file',k) obtains  subsequent  MATLAB  input  from  an
  1132.       external  file.  The printing of input is controlled by the
  1133.       optional parameter k .
  1134.       If k = 1 , the input is echoed.
  1135.       If k = 2 , the MATLAB prompt <> is printed.
  1136.       If k = 4 , MATLAB pauses before each prompt and waits for a
  1137.       null line to continue.
  1138.       If k = 0 , there is no echo, prompt or pause.  This is  the
  1139.       default if the exec command is followed by a semicolon.
  1140.       If k = 7 , there will be echos, prompts and pauses. This is
  1141.       useful for demonstrations on video terminals.
  1142.       If k = 3 , there will be echos and prompts, but no  pauses.
  1143.       This is the the default if the exec command is not followed
  1144.       by a semicolon.
  1145.       EXEC(0) causes subsequent input to  be  obtained  from  the
  1146.       terminal. An end-of-file has the same effect.
  1147.       EXEC's may be nested, i.e. the text in the file may contain
  1148.       EXEC of another file.  EXEC's may also be driven by FOR and
  1149.       WHILE loops.
  1150.  
  1151.  
  1152.  
  1153.  
  1154.  
  1155.  
  1156.  
  1157.  
  1158.  
  1159. MATLAB, page 52
  1160.  
  1161.  
  1162.  
  1163. EXIT  Causes termination of a FOR or WHILE loop.
  1164.       If not in a loop, terminates execution of MATLAB.
  1165.  
  1166. EXP   EXP(X)  is the exponential of  X ,  e  to the X .  See  FUN
  1167.       .
  1168.  
  1169. EYE   Identity matrix.  EYE(N) is the N  by  N  identity  matrix.
  1170.       EYE(M,N)   is an M by N matrix with 1's on the diagonal and
  1171.       zeros elsewhere.  EYE(A)  is the same size  as   A  .   EYE
  1172.       with  no  arguments is an identity matrix of whatever order
  1173.       is appropriate in the context.   For  example,  A  +  3*EYE
  1174.       adds  3  to each diagonal element of  A .
  1175.  
  1176. FILE  The EXEC, SAVE, LOAD,  PRINT  and  DIARY  functions  access
  1177.       files.   The  'file'  parameter  takes  different forms for
  1178.       different operating systems.  On most systems,  'file'  may
  1179.       be a string of up to 32 characters in quotes.  For example,
  1180.       SAVE('A') or EXEC('matlab/demo.exec') .  The string will be
  1181.       used as the name of a file in the local operating system.
  1182.       On all systems, 'file' may be a positive integer   k   less
  1183.       than  10  which  will  be  used  as  a FORTRAN logical unit
  1184.       number. Some systems then automatically access a file  with
  1185.       a  name  like  FORT.k  or FORk.DAT. Other systems require a
  1186.       file with a name like FT0kF001 to be assigned  to  unit   k
  1187.       before  MATLAB  is  executed. Check your local installation
  1188.       for details.
  1189.  
  1190. FLOPS Count of floating point operations.
  1191.       FLOPS  is  a  permanently  defined  row  vector  with   two
  1192.       elements.    FLOPS(1)  is  the  number  of  floating  point
  1193.       operations counted during the previous statement.  FLOPS(2)
  1194.       is  a  cumulative total.  FLOPS can be used in the same way
  1195.       as any other vector.  FLOPS(2) = 0  resets  the  cumulative
  1196.       total.   In  addition,  FLOPS(1) will be printed whenever a
  1197.       statement is terminated by an extra comma.  For example,
  1198.       X = INV(A);,
  1199.       or
  1200.       COND(A),   (as the last statement on the line).
  1201.       HELP FLPS gives more details.
  1202.  
  1203. FLPS  More detail on FLOPS.
  1204.       It is not feasible to count absolutely all  floating  point
  1205.       operations,  but  most  of  the important ones are counted.
  1206.       Each multiply and add in a real vector operation such as  a
  1207.       dot  product  or  a 'saxpy' counts one flop.  Each multiply
  1208.       and add in a complex vector  operation  counts  two  flops.
  1209.       Other additions, subtractions and multiplications count one
  1210.       flop each if the result is real and two flops if it is not.
  1211.       Real  divisions  count one and complex divisions count two.
  1212.       Elementary functions count one if real and two if  complex.
  1213.       Some examples.  If A and B are real N by N matrices, then
  1214.       A + B  counts N**2 flops,
  1215.       A*B    counts N**3 flops,
  1216.  
  1217.  
  1218.  
  1219.  
  1220.  
  1221.  
  1222.  
  1223.  
  1224.  
  1225. MATLAB, page 53
  1226.  
  1227.  
  1228.  
  1229.       A**100 counts 99*N**3 flops,
  1230.       LU(A)  counts roughly (1/3)*N**3 flops.
  1231.  
  1232. FOR   Repeat statements a specific number of times.
  1233.       FOR variable = expr, statement, ..., statement, END
  1234.       The END at the end of a line may  be  omitted.   The  comma
  1235.       before  the  END  may  also be omitted.  The columns of the
  1236.       expression are stored one at a time  in  the  variable  and
  1237.       then the following statements, up to the END, are executed.
  1238.       The expression is often of the form X:Y, in which case  its
  1239.       columns  are  simply  scalars.  Some examples (assume N has
  1240.       already been assigned a value).
  1241.       FOR I = 1:N, FOR J = 1:N, A(I,J) = 1/(I+J-1);
  1242.       FOR J = 2:N-1, A(J,J) = J; END; A
  1243.       FOR S = 1.0: -0.1: 0.0, ...  steps S with increments of -0.1 .
  1244.       FOR E = EYE(N), ...   sets  E  to the unit N-vectors.
  1245.       FOR V = A, ...   has the same effect as
  1246.       FOR J = 1:N, V = A(:,J); ...  except J is also set here.
  1247.  
  1248. FUN   For matrix arguments  X , the  functions  SIN,  COS,  ATAN,
  1249.       SQRT,  LOG,  EXP and X**p are computed using eigenvalues  D
  1250.       and eigenvectors  V .  If  <V,D> =  EIG(X)   then   f(X)  =
  1251.       V*f(D)/V  .   This method may give inaccurate results if  V
  1252.       is badly conditioned.  Some idea of  the  accuracy  can  be
  1253.       obtained by comparing  X**1  with  X .
  1254.       For vector arguments,  the  function  is  applied  to  each
  1255.       component.
  1256.  
  1257. HESS  Hessenberg form.  The Hessenberg form of a matrix  is  zero
  1258.       below the first subdiagonal.  If the matrix is symmetric or
  1259.       Hermitian,  the  form  is  tridiagonal.   <P,H>  =  HESS(A)
  1260.       produces  a  unitary  matrix P and a Hessenberg matrix H so
  1261.       that A = P*H*P'.  By itself, HESS(A) returns H.
  1262.  
  1263. HILB  Inverse Hilbert matrix.  HILB(N)  is the inverse of  the  N
  1264.       by  N   matrix  with elements  1/(i+j-1), which is a famous
  1265.       example of a badly conditioned matrix.  The result is exact
  1266.       for  N  less than about 15, depending upon the computer.
  1267.  
  1268. IF    Conditionally execute statements.  Simple form...
  1269.       IF expression rop expression, statements
  1270.       where rop is =, <, >, <=, >=, or  <>  (not  equal)  .   The
  1271.       statements  are  executed  once if the indicated comparison
  1272.       between the real parts of the first components of  the  two
  1273.       expressions  is true, otherwise the statements are skipped.
  1274.       Example.
  1275.       IF ABS(I-J) = 1, A(I,J) = -1;
  1276.       More complicated forms use END in the same way it  is  used
  1277.       with FOR and WHILE and use ELSE as an abbreviation for END,
  1278.       IF expression not rop expression .  Example
  1279.       FOR I = 1:N, FOR J = 1:N, ...
  1280.          IF I = J, A(I,J) = 2; ELSE IF ABS(I-J) = 1, A(I,J) = -1; ...
  1281.          ELSE A(I,J) = 0;
  1282.  
  1283.  
  1284.  
  1285.  
  1286.  
  1287.  
  1288.  
  1289.  
  1290.  
  1291. MATLAB, page 54
  1292.  
  1293.  
  1294.  
  1295.       An easier way to accomplish the same thing is
  1296.       A = 2*EYE(N);
  1297.       FOR I = 1:N-1, A(I,I+1) = -1; A(I+1,I) = -1;
  1298.  
  1299. IMAG  IMAG(X)  is the imaginary part of  X .
  1300.  
  1301. INV   INV(X)  is the inverse of the square matrix  X .  A warning
  1302.       message  is  printed  if   X   is  badly  scaled  or nearly
  1303.       singular.
  1304.  
  1305. KRON  KRON(X,Y) is the Kronecker tensor product of X and Y  .  It
  1306.       is  also  denoted by X .*. Y . The result is a large matrix
  1307.       formed by taking all possible products between the elements
  1308.       of  X  and  those  of Y . For example, if X is 2 by 3, then
  1309.       X .*. Y is
  1310.  
  1311.             < x(1,1)*Y  x(1,2)*Y  x(1,3)*Y
  1312.               x(2,1)*Y  x(2,2)*Y  x(2,3)*Y >
  1313.  
  1314.       The five-point discrete Laplacian for an n-by-n grid can be
  1315.       generated by
  1316.  
  1317.             T = diag(ones(n-1,1),1);  T = T + T';  I = EYE(T);
  1318.             A = T.*.I + I.*.T - 4*EYE;
  1319.  
  1320.       Just  in  case  they  might  be  useful,  MATLAB   includes
  1321.       constructions called Kronecker tensor quotients, denoted by
  1322.       X ./. Y and X .\. Y .  They are obtained by  replacing  the
  1323.       elementwise multiplications in X .*. Y with divisions.
  1324.  
  1325. LINES An internal count is kept of the number of lines of  output
  1326.       since  the  last  input.   Whenever this count approaches a
  1327.       limit, the  user  is  asked  whether  or  not  to  suppress
  1328.       printing  until the next input.  Initially the limit is 25.
  1329.       LINES(N) resets the limit to N .
  1330.  
  1331. LOAD  LOAD('file') retrieves all the variables from  the  file  .
  1332.       See  FILE  and  SAVE for more details.  To prepare your own
  1333.       file for LOADing, change the READs to WRITEs  in  the  code
  1334.       given under SAVE.
  1335.  
  1336. LOG   LOG(X)  is the  natural  logarithm  of   X  .   See  FUN  .
  1337.       Complex results are produced if  X  is not positive, or has
  1338.       nonpositive eigenvalues.
  1339.  
  1340. LONG  Determine output format.   All  computations  are  done  in
  1341.       complex arithmetic and double precision if it is available.
  1342.       SHORT and  LONG  merely  switch  between  different  output
  1343.       formats.
  1344.       SHORT    Scaled fixed point format with about 5 digits.
  1345.       LONG     Scaled fixed point format with about 15 digits.
  1346.       SHORT E  Floating point format with about 5 digits.
  1347.       LONG E   Floating point format with about 15 digits.
  1348.  
  1349.  
  1350.  
  1351.  
  1352.  
  1353.  
  1354.  
  1355.  
  1356.  
  1357. MATLAB, page 55
  1358.  
  1359.  
  1360.  
  1361.       LONG Z   System dependent format, often hexadecimal.
  1362.  
  1363. LU    Factors from Gaussian elimination.  <L,U> = LU(X)  stores a
  1364.       upper triangular matrix in  U  and a 'psychologically lower
  1365.       triangular matrix', i.e. a product of lower triangular  and
  1366.       permutation matrices, in L , so that  X = L*U .  By itself,
  1367.       LU(X) returns the output from CGEFA .
  1368.  
  1369. MACRO The macro facility involves text and inward pointing  angle
  1370.       brackets.  If  STRING  is  the  source  text for any MATLAB
  1371.       expression or statement, then
  1372.             t = 'STRING';
  1373.       encodes the text as a vector of integers  and  stores  that
  1374.       vector in  t .  DISP(t) will print the text and
  1375.             >t<
  1376.       causes the text to be interpreted, either as a statement or
  1377.       as a factor in an expression.  For example
  1378.             t = '1/(i+j-1)';
  1379.             disp(t)
  1380.             for i = 1:n, for j = 1:n, a(i,j) = >t<;
  1381.       generates the Hilbert matrix of order n.
  1382.       Another example showing indexed text,
  1383.             S = <'x = 3            '
  1384.                  'y = 4            '
  1385.                  'z = sqrt(x*x+y*y)'>
  1386.             for k = 1:3, >S(k,:)<
  1387.       It is necessary that the strings making up  the  "rows"  of
  1388.       the "matrix"  S  have the same lengths.
  1389.  
  1390. MAGIC Magic square.  MAGIC(N) is an N  by  N  matrix  constructed
  1391.       from  the integers 1 through N**2 with equal row and column
  1392.       sums.
  1393.  
  1394. NORM  For matrices..
  1395.       NORM(X)  is the largest singular value of  X .
  1396.       NORM(X,1)  is the 1-norm of  X .
  1397.       NORM(X,2)  is the same as NORM(X) .
  1398.       NORM(X,'INF')  is the infinity norm of  X .
  1399.       NORM(X,'FRO')  is the F-norm, i.e.  SQRT(SUM(DIAG(X'*X))) .
  1400.       For vectors..
  1401.       NORM(V,P) = (SUM(V(I)**P))**(1/P) .
  1402.       NORM(V) = NORM(V,2) .
  1403.       NORM(V,'INF') = MAX(ABS(V(I))) .
  1404.  
  1405. ONES  All ones.  ONES(N)  is an N by N matrix of ones.  ONES(M,N)
  1406.       is an M by N matrix of ones .  ONES(A)  is the same size as
  1407.       A  and all ones .
  1408.  
  1409. ORTH  Orthogonalization.   Q  =  ORTH(X)   is   a   matrix   with
  1410.       orthonormal  columns,  i.e. Q'*Q = EYE, which span the same
  1411.       space as the columns of  X .
  1412.  
  1413. PINV  Pseudoinverse.  X = PINV(A) produces a matrix   X   of  the
  1414.  
  1415.  
  1416.  
  1417.  
  1418.  
  1419.  
  1420.  
  1421.  
  1422.  
  1423. MATLAB, page 56
  1424.  
  1425.  
  1426.  
  1427.       same  dimensions as  A' so that  A*X*A = A , X*A*X = X  and
  1428.       AX  and  XA  are Hermitian .  The computation is  based  on
  1429.       SVD(A)  and  any  singular values less than a tolerance are
  1430.       treated   as    zero.     The    default    tolerance    is
  1431.       NORM(SIZE(A),'inf')*NORM(A)*EPS.   This  tolerance  may  be
  1432.       overridden with X = PINV(A,tol).  See RANK.
  1433.  
  1434. PLOT  PLOT(X,Y) produces a plot of  the  elements  of  Y  against
  1435.       those  of X . PLOT(Y) is the same as PLOT(1:n,Y) where n is
  1436.       the  number  of   elements   in   Y   .    PLOT(X,Y,P)   or
  1437.       PLOT(X,Y,p1,...,pk)  passes the optional parameter vector P
  1438.       or scalars p1 through pk to the plot routine.  The  default
  1439.       plot  routine  is a crude printer-plot. It is hoped that an
  1440.       interface to local graphics equipment can be provided.
  1441.       An interesting example is
  1442.             t = 0:50;
  1443.             PLOT( t.*cos(t), t.*sin(t) )
  1444.  
  1445. POLY  Characteristic polynomial.
  1446.       If  A  is an N by N matrix, POLY(A) is a column vector with
  1447.       N+1   elements   which   are   the   coefficients   of  the
  1448.       characteristic polynomial,  DET(lambda*EYE - A) .
  1449.       If V is a vector, POLY(V) is a vector  whose  elements  are
  1450.       the  coefficients  of  the  polynomial  whose roots are the
  1451.       elements of V .  For vectors, ROOTS and  POLY  are  inverse
  1452.       functions  of  each  other,  up  to  ordering, scaling, and
  1453.       roundoff error.
  1454.       ROOTS(POLY(1:20)) generates Wilkinson's famous example.
  1455.  
  1456. PRINT PRINT('file',X) prints X on  the  file  using  the  current
  1457.       format determined by SHORT, LONG Z, etc.  See FILE.
  1458.  
  1459. PROD  PROD(X)  is the product of all the elements of  X .
  1460.  
  1461. QR    Orthogonal-triangular decomposition.
  1462.       <Q,R> = QR(X)  produces an upper triangular  matrix   R  of
  1463.       the  same dimension as  X  and a unitary matrix  Q  so that
  1464.       X = Q*R .
  1465.       <Q,R,E> = QR(X)  produces a  permutation  matrix   E  ,  an
  1466.       upper  triangular  R  with decreasing diagonal elements and
  1467.       a unitary  Q  so that  X*E = Q*R .
  1468.       By itself, QR(X) returns the output of CQRDC .  TRIU(QR(X))
  1469.       is R .
  1470.  
  1471. RAND  Random numbers and matrices.  RAND(N)  is an N by N  matrix
  1472.       with  random  entries.  RAND(M,N)  is an M by N matrix with
  1473.       random entries.  RAND(A)  is the same size as   A  .   RAND
  1474.       with no arguments is a scalar whose value changes each time
  1475.       it is referenced.
  1476.       Ordinarily,  random numbers are  uniformly  distributed  in
  1477.       the  interval  (0.0,1.0)  .   RAND('NORMAL')  switches to a
  1478.       normal distribution  with  mean  0.0  and  variance  1.0  .
  1479.       RAND('UNIFORM')  switches back to the uniform distribution.
  1480.  
  1481.  
  1482.  
  1483.  
  1484.  
  1485.  
  1486.  
  1487.  
  1488.  
  1489. MATLAB, page 57
  1490.  
  1491.  
  1492.  
  1493.       RAND('SEED') returns the current value of the seed for  the
  1494.       generator.    RAND('SEED',n)   sets   the   seed   to  n  .
  1495.       RAND('SEED',0) resets the seed to 0, its value when  MATLAB
  1496.       is first entered.
  1497.  
  1498. RANK  Rank.  K = RANK(X) is the number of singular values  of   X
  1499.       that are larger than NORM(SIZE(X),'inf')*NORM(X)*EPS.
  1500.       K = RANK(X,tol) is the number of singular values of  X that
  1501.       are larger than tol .
  1502.  
  1503. RCOND RCOND(X)   is  an  estimate  for  the  reciprocal  of   the
  1504.       condition  of   X   in  the  1-norm obtained by the LINPACK
  1505.       condition estimator.  If  X  is well conditioned,  RCOND(X)
  1506.       is  near  1.0  .   If  X  is badly conditioned, RCOND(X) is
  1507.       near 0.0 .
  1508.       <R, Z> = RCOND(A) sets  R  to RCOND(A) and also produces  a
  1509.       vector  Z so that
  1510.                  NORM(A*Z,1) = R*NORM(A,1)*NORM(Z,1)
  1511.       So, if RCOND(A) is small, then  Z  is an  approximate  null
  1512.       vector.
  1513.  
  1514. RAT   An experimental  function  which  attempts  to  remove  the
  1515.       roundoff   error  from  results  that  should  be  "simple"
  1516.       rational numbers.
  1517.       RAT(X) approximates each  element  of   X  by  a  continued
  1518.       fraction of the form
  1519.  
  1520.                 a/b = d1 + 1/(d2 + 1/(d3 + ... + 1/dk))
  1521.  
  1522.       with k <= len, integer di and abs(di) <= max .  The default
  1523.       values of the parameters are len = 5 and max = 100.
  1524.       RAT(len,max) changes the default values.  Increasing either
  1525.       len or max increases the number of possible fractions.
  1526.       <A,B> = RAT(X) produces integer matrices A and B so that
  1527.  
  1528.                 A ./ B  =  RAT(X)
  1529.  
  1530.       Some examples:
  1531.  
  1532.             long
  1533.             T = hilb(6), X = inv(T)
  1534.             <A,B> = rat(X)
  1535.             H = A ./ B, S = inv(H)
  1536.  
  1537.             short e
  1538.             d = 1:8,  e = ones(d),  A = abs(d'*e - e'*d)
  1539.             X = inv(A)
  1540.             rat(X)
  1541.             display(ans)
  1542.  
  1543.  
  1544. REAL  REAL(X)  is the real part of  X .
  1545.  
  1546.  
  1547.  
  1548.  
  1549.  
  1550.  
  1551.  
  1552.  
  1553.  
  1554.  
  1555. MATLAB, page 58
  1556.  
  1557.  
  1558.  
  1559. RETURN  From the terminal, causes return to the operating  system
  1560.       or  other  program  which  invoked  MATLAB.  From inside an
  1561.       EXEC, causes  return  to  the  invoking  EXEC,  or  to  the
  1562.       terminal.
  1563.  
  1564. RREF  RREF(A) is the reduced row echelon form of the  rectangular
  1565.       matrix.  RREF(A,B) is the same as RREF(<A,B>) .
  1566.  
  1567. ROOTS Find polynomial roots.  ROOTS(C)  computes the roots of the
  1568.       polynomial  whose  coefficients  are  the  elements  of the
  1569.       vector  C .  If  C  has  N+1  components, the polynomial is
  1570.       C(1)*X**N + ... + C(N)*X + C(N+1) .  See POLY.
  1571.  
  1572. ROUND ROUND(X)  rounds  the  elements  of   X   to  the   nearest
  1573.       integers.
  1574.  
  1575. SAVE  SAVE('file') stores all the current variables in a file.
  1576.       SAVE('file',X) saves only X .  See FILE .
  1577.       The variables may be retrieved later by LOAD('file') or  by
  1578.       your  own program using the following code for each matrix.
  1579.       The lines involving XIMAG may be eliminated  if  everything
  1580.       is known to be real.
  1581.  
  1582.             attach lunit to 'file'
  1583.             REAL or DOUBLE PRECISION XREAL(MMAX,NMAX)
  1584.             REAL or DOUBLE PRECISION XIMAG(MMAX,NMAX)
  1585.             READ(lunit,101) ID,M,N,IMG
  1586.             DO 10 J = 1, N
  1587.                READ(lunit,102) (XREAL(I,J), I=1,M)
  1588.                IF (IMG .NE. 0) READ(lunit,102) (XIMAG(I,J),I=1,M)
  1589.          10 CONTINUE
  1590.  
  1591.       The formats used are system dependent.  The  following  are
  1592.       typical.     See    SUBROUTINE   SAVLOD   in   your   local
  1593.       implementation of MATLAB.
  1594.  
  1595.         101 FORMAT(4A1,3I4)
  1596.         102 FORMAT(4Z18)
  1597.         102 FORMAT(4O20)
  1598.         102 FORMAT(4D25.18)
  1599.  
  1600. SCHUR Schur decomposition.  <U,T> = SCHUR(X)  produces  an  upper
  1601.       triangular  matrix   T , with the eigenvalues of  X  on the
  1602.       diagonal, and a unitary matrix  U so that  X =  U*T*U'  and
  1603.       U'*U = EYE .  By itself, SCHUR(X) returns  T .
  1604.  
  1605. SHORT See LONG .
  1606.  
  1607. SEMI  Semicolons at the end of  lines  will  cause,  rather  than
  1608.       suppress,  printing.   A  second  SEMI restores the initial
  1609.       interpretation.
  1610.  
  1611. SIN   SIN(X)  is the sine of  X .  See FUN .
  1612.  
  1613.  
  1614.  
  1615.  
  1616.  
  1617.  
  1618.  
  1619.  
  1620.  
  1621. MATLAB, page 59
  1622.  
  1623.  
  1624.  
  1625. SIZE  If X is an M by N matrix, then SIZE(X) is <M, N> .
  1626.       Can also be used with a multiple assignment,
  1627.             <M, N> = SIZE(X) .
  1628.  
  1629. SQRT  SQRT(X)  is the square root of  X .   See  FUN  .   Complex
  1630.       results  are  produced  if   X   is  not  positive,  or has
  1631.       nonpositive eigenvalues.
  1632.  
  1633. STOP  Use EXIT instead.
  1634.  
  1635. SUM   SUM(X)   is  the  sum  of  all  the  elements   of    X   .
  1636.       SUM(DIAG(X))  is the trace of  X .
  1637.  
  1638. SVD   Singular value decomposition.  <U,S,V> = SVD(X)  produces a
  1639.       diagonal  matrix  S , of the same dimension as  X  and with
  1640.       nonnegative diagonal  elements  in  decreasing  order,  and
  1641.       unitary matrices  U  and  V  so that  X = U*S*V' .
  1642.       By itself, SVD(X) returns a vector containing the  singular
  1643.       values.
  1644.       <U,S,V>   =   SVD(X,0)   produces   the   "economy    size"
  1645.       decomposition.   If  X  is m by n with m > n, then only the
  1646.       first n columns of U are computed and S is n by n .
  1647.  
  1648. TRIL  Lower triangle.  TRIL(X) is the lower triangular part of X.
  1649.       TRIL(X,K) is the elements on and below the K-th diagonal of
  1650.       X.  K = 0 is the main diagonal, K > 0  is  above  the  main
  1651.       diagonal and K < 0 is below the main diagonal.
  1652.  
  1653. TRIU  Upper triangle.  TRIU(X) is the upper triangular part of X.
  1654.       TRIU(X,K) is the elements on and above the K-th diagonal of
  1655.       X.  K = 0 is the main diagonal, K > 0  is  above  the  main
  1656.       diagonal and K < 0 is below the main diagonal.
  1657.  
  1658. USER  Allows personal  Fortran  subroutines  to  be  linked  into
  1659.       MATLAB .  The subroutine should have the heading
  1660.  
  1661.                SUBROUTINE USER(A,M,N,S,T)
  1662.                REAL or DOUBLE PRECISION A(M,N),S,T
  1663.  
  1664.       The MATLAB statement  Y = USER(X,s,t)  results in a call to
  1665.       the  subroutine with a copy of the matrix  X  stored in the
  1666.       argument  A , its column and row dimensions in  M  and  N ,
  1667.       and  the scalar parameters  s  and  t  stored in  S  and  T
  1668.       . If  s and t  are omitted, they are set to  0.0  .   After
  1669.       the  return,   A  is stored in  Y .  The dimensions  M  and
  1670.       N  may be reset within the subroutine.  The statement  Y  =
  1671.       USER(K)  results in a call with M = 1, N = 1  and  A(1,1) =
  1672.       FLOAT(K) .  After the subroutine has been written, it  must
  1673.       be compiled and linked to the MATLAB object code within the
  1674.       local operating system.
  1675.  
  1676. WHAT  Lists commands and functions currently available.
  1677.  
  1678.  
  1679.  
  1680.  
  1681.  
  1682.  
  1683.  
  1684.  
  1685.  
  1686.  
  1687. MATLAB, page 60
  1688.  
  1689.  
  1690.  
  1691. WHILE Repeat statements an indefinite number of times.
  1692.       WHILE expr rop expr, statement, ..., statement, END
  1693.       where rop is =, <, >, <=, >=, or <> (not equal) .  The  END
  1694.       at  the end of a line may be omitted.  The comma before the
  1695.       END may also be omitted.  The commas  may  be  replaced  by
  1696.       semicolons   to   avoid   printing.    The  statements  are
  1697.       repeatedly executed as long  as  the  indicated  comparison
  1698.       between  the  real parts of the first components of the two
  1699.       expressions is true.   Example  (assume  a  matrix   A   is
  1700.       already defined).
  1701.       E = 0*A; F = E + EYE; N = 1;
  1702.       WHILE NORM(E+F-E,1) > 0, E = E + F; F = A*F/N; N = N + 1;
  1703.       E
  1704.  
  1705. WHO   Lists current variables.
  1706.  
  1707. WHY   Provides succinct answers to any questions.
  1708.  
  1709. //
  1710.  
  1711.  
  1712.  
  1713.  
  1714.  
  1715.  
  1716.  
  1717.  
  1718.  
  1719.  
  1720.  
  1721.  
  1722.  
  1723.  
  1724.  
  1725.  
  1726.  
  1727.  
  1728.  
  1729.  
  1730.  
  1731.  
  1732.  
  1733.  
  1734.  
  1735.  
  1736.  
  1737.  
  1738.  
  1739.  
  1740.  
  1741.  
  1742.  
  1743.  
  1744.  
  1745.  
  1746.  
  1747.  
  1748.  
  1749.  
  1750.  
  1751. SHAR_EOF
  1752. #    End of shell archive
  1753. exit 0
  1754. -- 
  1755. Bob Page, U of Lowell CS Dept.  page@swan.ulowell.edu  ulowell!page
  1756. Have five nice days.
  1757.